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The coefficient of diffusion of hydrogen in crystalline silicon is calculated using tight-binding molec- 
ular dynamics. Our results are in good quantitative agreement with an earlier study by Panzarini 
and Colombo [Phys. Rev. Lett. 73, 1636 (1994)]. However, while our calculations indicate that long 
jumps dominate over single hops at high temperatures, no abrupt change in the diffusion coefficient 
can be observed with decreasing temperature. The (classical) Arrhenius diffusion parameters, as a 
consequence, should extrapolate to low temperatures. 

PACS numbers: 66.30.Jt, 66.30. Dn 
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In spite of the tremendous efforts that have been en- 
gaged in determining the rate of diffusion of hydrogen in 
crystalline silicon, a concensus on the "true value" has 
not yet been reached. Of course, the diffusion constant 
varies strongly with temperature — not necessarily in a 
perfect Arrhenius manner — making a precise determina- 
tion of the diffusion parameters difficult. Other complica- 
tions arise from possible collective effects (as opposed to 
tracer diffusion), low-temperature quantum effects, im- 
purities and defects, etc. 

Experimental estimates afp-the diffusion constant are 
numerous and vary widelyuLl sometimes by two orders 
of magnitude at a given temperature, as can be appre- 
ciated from the open circles in Fig. |l](a). This is an un- 
pleasant state of affairs, since precise knowledge of this 
quantity is important for both practical and fundamental 
reasons: Because it forms complexes with a variety of de- 
fects, hydrogen affects deeply the optical and electronic 
properties of semiconductors. It is usually present as a re- 
sult of the fabrication process, but is often intentionally 
introduced in order to passivate defects. Being a light 
species, further, H diffuses readily, inducing additional 
defects along its way, thus affecting the transport prop- 
erties of the material to an extent which is determined 
by its relative concentration. It is therefore important 
to understand diffusion at the atomic level in order to 
gain better control on the properties of semiconductors 
and, in view of the simplicity of the structure of the host 
material, it is of fundamental importance to be able to 
understand this prototypical system. 

The diffusion coefficient can be estimated, at suffi- 
ciently high temperature, using the now well-established 
molecular-dynamics method; a proper model for the in- 
teratomic potentials is then needed. Empirical poten- 
tials lack the transferability and predictive power of first- 
principles methods. The latter, however, are subject to 
limitations in size and time, which makes them unpracti- 



cal for the long simulations required for a proper (statis- 
tically meaningful) estimate of the rate of diffusion. In 
an early application of the Car-Parrinello met hod, El Buda 
et alB calculated the coefficient of diffusion of H + in Si at 
three temperatures in the range 1200-1950 K, covering a 
minuscule maximum observation time of 4 ps. 

The semi-empirical, semi-quantum 

tight-binding molecular-dynamics (TBMD) scheme, orig- 
inally proposed by Khan and BroughtonE] and Goodwin, 
Skinner and Pettifor (GSP)J12l provideSj-good accuracy 
at a very reasonable computational cost.Ej Here, the at- 
tractive part of the atom-atom interactions is described 
quantum-mechanically using (parametrized) overlap (or 
hopping) integrals, while the repulsive part is fitted from 
known properties of the system, e.g., binding energy vs 
distance. The forces are then derived from the Hellman- 
Feynman theorem. TBMD models for_Si:H were pro- 
posed by Panzarini and Colombo (PCUp Boucher and 
DeLeo (BDL),EJ and Kim, Lee and Leejia all three mod- 
els are based on the GSP model for Si-Si interactions, 
and use comparable fitting schemes to take into account 
Si-H and H-H interactions. 

The problem of hydrogen diffusion in silicon—was ad- 
dressed using TBMD by Paazarini and ColombcEj as well 
as by Boucher and DeLeo;EJ both considered a single H 
atom in a 64-atom c-Si supercell. While extending sig- 
nificantly the range of temperatures that .were covered 
by the ah initio simulations of Buda et alB (1050-2000 
K for BDL, 800-1800 K for PC), the timescales covered 
by these MD calculations remain short (42 and 300 ps, 
respectively): indeed, for a diffusion constant of I0~ 6 
cm 2 /s — as found at about 800 K — a quick calcula- 
tion indicates that the average size of the region visited 
by the diffusing particle over 300 ps would be about 4 
A, corresponding roughly to the second-neighbour dis- 
tance in c-Si. Further, the two calculations exhibit some 
disagreement which might be inherent to the models or, 
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more likely, to the statistical quality of the MD data. 

In this short note, we revisit the problem using, again, 
TBMD (PC version), but with much longer timescales: 
our simulations ran during a formidable 7 nanoseconds at 
the lowest temperature we could decently examine — 700 
K, which is 100 K below the lowest temperature looked at 
by PC. Our calculations generally confirm PC's results, 
in particular the discrepancies with experiment observed 
at the lowest temperatures. However, while our calcula- 
tions indicate that long jumps dominate over single hops 
at high temperatures, the diffusion coefficient exhibits no 
abrupt change with temperature. We are led to conclude 
that the Arrhenius diffusion parameters should extrapo- 
late to low temperatures as far as the classical part of the 
motion is concerned. 

As previewed earlier, Fig. |l](a) presents the results of 
several measurements of the diffusion constant, plotted a 
la Arrhenius. r-Also indicated are the ab initio MD data 
of Buda et aZ.Jj they are found to be in reasonable agree- 
ment with the high-temperatuxe experimental points of 
Van Wieringen and WarmoltzEl, fitted to the Arrhenius 
law D(T) = D Q exp(-E A /k B T), with D Q = 9.41 x 10~ 3 
cm 2 /s and E A — 0.48 eV. When extended to low tem- 
peratures [dotted line in Fig. H(a)], one clearly sees the 
deviations from the Arrhenius behaviour; it should be 
said, however, that there is no "guarantee" that diffu- 
sion should be Arrhenius over the whole range of tem- 
peratures. 

The TBMD data of PC and BDL are also plotted in 
Fig. 0(a), and more legibly in Fig. 0(b). The agreement 
with experiment is clearly excellent at high temperature 
— certainly within the errors that can be associated with 
both measurements and calculations. The data of BDL 
are found to be extremely well fitted by the Arrhenius 
law with D = 6.91 x I0~ 3 cm 2 /s and E A = 0.45 eV 
all the way down to 1050 K, in striking agreement with 
experiment, as can be judged by the close similarity be- 
tween the prefactors and energy barriers. PC, in contrast, 
observe significant deviations from Arrhenius already at 
1200 K, a problem that can possibly be attributed to the 
statistical quality of their data. 

In Fig. ^ we give, as an example, the calculated time- 
dependence of the mean-square displacement for our 
lowest-temperature run, viz. 700 K. The simulation at 
this temperature ran for a total of 7 ns and the mean- 
square displacement is calculated for a maximum cor- 
relation time of 1 ns in order to minimize statistical 
uncertainties. We also evaluated the mean-square dis- 
placement using only the first half of the run, then 
only the second half so as to have a feeling for the 
"error bar" of our estimate of the diffusion coefficient, 
D = Hindoo r 2 (t)/6t. The three different calculations 
are indicated in Fig. || by full, dashed, and dotted lines, 
respectively. The corresponding coefficients of diffusion, 
which we discuss next, are presented in Fig. [IJ for each 
temperature we investigated (700, 800, 900, 1000, 1200, 
and 1500 K), three data points are thus given. 

As can be appreciated from Fig. 0(b), our estimates of 



the diffusion coefficient generally agree with the values 
reported by PC, but clearly with much improved statis- 
tical quality: the data points exhibit very little scatter, 
falling along a single, rather well-defined straight line. 
It is of course hazardous to assume that the diffusion 
process is perfectly Arrhenius and that no "exotic" diffu- 
sion mechanisms (i.e., other than single hops) are taking 
place. In order to "guide the eye" , however, we fitted our 
data to the law D = 8.9 x I0" 3 exp(0.58eV//c B T) cm 2 /s; 
this is displayed as the dashed line in Fig. |[ Our calcu- 
lations, clearly, show no evidence of a sizeable change in 
the diffusive behaviour as a function of temperature, as 
can possibly be (and was) inferred from the data of PC. 

Exotic mechanisms — in the present case long jumps 
— are in fact present; this has already been noted by PC, 
who also found long jumps to be quenched in as tem- 
perature drops; this observation was based on a visual 
examination of the mean-square displacements. A phys- 
ically more appropriate characterization of single-atom 
motion is provided by the self part of the van Hove cor- 
relation function (see, e.g., Ref. |f^), G s (r,t), which gives 
the probability of finding a particle at r at time t given 
that it was at the origin at t = 0. Averaging out over 
angular space, the function of interest is 47rr 2 G s (r, t). 

The van Hove self-correlation function for the H atom 
in c-Si at 700 K is displayed in Fig. || Here, 47rr 2 G s (r, t) 
is plotted as a function of distance for four different times. 
At short time (but somewhat longer than the typical vi- 
brational period), single hops dominate the motion. As 
time increases, the particle is (of course) on average fur- 
ther away and, evidently, the diffusive motion proceeds 
via a sequence of single jumps — to a reasonable approx- 
imation. 

Because diffusion is activated, it is more meaningful, 
in order to compare the modes of diffusion at different 
temperatures, to examine the behaviour of 4irr 2 G s (r, t) 
at fixed mean-square displacement. Thus, it is possible 
in this way to determine the type of motion that leads 
a particle a certain distance away from its original posi- 
tion. Here we chose (somewhat arbitrarily) a fixed mean- 
square displacement of 50 A 2 . From the r 2 (t) curves, it 
is easy to determine the time at which, on average, and 
for each temperature, the particle will be at the required 
fixed position. 

The results of this analysis are presented in Fig. ||. 
Again, here, we find that the motion consists of single 
hops at low temperatures, but acquires a "long-jump" 
character as temperature increases. At the highest tem- 
peratures, in fact, the distribution is nearly continuous, 
and single jumps become almost undetectable. There is 
some evidence, in these plots, of a discontinuous change 
in the manner that diffusion proceeds as a function of 
temperature: upon going from 800 to 900 K, nearest- 
neighbour jumps are found to become negligible, to the 
advantage of second- and further- neighbour hops. But 
then no sign of such a change is apparent in the diffu- 
sion coefficient (cf. Fig. l|), which perhaps lacks the sen- 
sitivity needed to reveal such details. Thus, we are led 
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to conclude that the low-temperature (T < about 900 
K) diffusion coefficient is indeed Arrhenius-like, whereas 
high-temperature diffusion proceeds via a complicated se- 
quence of jumps whose energy barriers combine into a 
single activation energy that is close to that for single 
hops. 

In view of this, there are no reasons to believe that this 
Arrhenius law should not extrapolate to very low temper- 
atures, if the system were classical. This is clearly not the 
case here and it is of course expected that diffusion will 
be enhanced by quantum contributions, even more so as 
temperature decreases. We have not included quantum 
corrections in the present study: this is a difficult prob- 
lem and it is not clear that their role is significant at the 
temperatures considered here. The small discrepancies 
between TBMD and experimental data at high temper- 
atures might be related to quantum effects, but are per- 
haps more likely due to some limitations of the model, 
i.e., are not significant. Likewise, Fig. |l| seems to indicate 
some partial "recovery" of the experimental data at very 
low temperature. This might be a manifestation, again, 
of quantum contributions. Clearly a more consistent set 
of experimental data is needed in order to have a proper 
perspective on the problem. 
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FIG. 1. (a) and (b) Coefficient of diffusion of hydrogen 
in crystalline silicon. Panel (b) zooms in on the high tem- 
perature region. The present data are indicated by the black 
squares; the dashed line is an Arrhenius fit to them. Other 
data are as follows: open circles — experimental data of Refs. 
0-0; dotted line — fit to the high-temperature data of Ref. 
[[land corresponding extrapolation to low temperatures; open 
squares — TBMD data of Panzarini and Colombo, Ref. |l3|; 

ren triangles: ab initio MD calculations of Buda el al. , Ref. 
solid line: fit to the TBMD results of Boucher and DeLeo, 
Ref. 0. 
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FIG. 2. Mean-square displacement of the hydrogen atom 
versus time at 700 K for a maximum correlation time of 1 ns 
out of a 7 ns run. The three lines correspond to averaging 
over the whole duration of the run (full line), the first half 
(dotted line) and the second half (dashed line). The differ- 
ences between the three curves give an idea of the error bar 
on the diffusion constant. 




FIG. 3. Van Hove self-correlation function at 700 K as a 
function of distance for four different correlation times: 1 ps 
(long dashes), 50 ps (short dashes), 200 ps (dots), and 1 ns 
(full line). 
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FIG. 4. Van Hove self-correlation function at various tem- 
peratures for fixed average mean-square displacement (50 A 2 ), 
as explained in the text. 
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